# R Script analyzes results from conjoint analysis
#
# author: Jonas Markgraf & Guillermo Rosas
# last updated: 2023-03-12
###############################################################################

rm(list = ls())

# load packages
# load libraries
packages <- c("cjoint", "repmis", "readxl", "tidyverse", "FindIt"
              , "stringr", "cregg", "sandwich", "lfe", "ggpubr")
new.packages <- packages[!(packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)
lapply(packages, require, character.only = TRUE)
rm(packages, new.packages)

# set wd
possibles <- c("~/Desktop/replication_files/") # set your wd here

set_valid_wd(possibles)

# source DiD functions
source("Conjoint_DiD-functions.R")

# load data
surveyData <- readRDS("cleanedConjointData.csv")

# Figure 2: AMCE : ========================================================

fig2 <- forcedChoice ~ mrtg.crdt + consum.crdt + welfare + unempl + tax
amceMain <- cj(surveyData, fig2, id = ~ResponseId
               , feature_labels = list(mrtg.crdt = "Housing Loan"
                                       , consum.crdt = "Credit Card Loan"
                                       , welfare = "Social Assistance"
                                       , unempl = "Unemployment Support"
                                       , tax = "Income Tax")
               , level_order = "descending")

#pdf(file = "graphPlots/AMCE_fullSample.pdf",
#    width = 7,
#    height = 3)
plot(amceMain, legend_pos = "none")
#dev.off()

# MAIN RESULTS (Table 3)  ========================================================

### getting numbers for Table 3 (compiled manually) ###

# Table 3, model 1 (full sample)
tab3m1q1 <- DifnDif3bis (a="tax", b="consum.crdt")
tab3m1q2 <- DifnDif3bis (a="tax", b="mrtg.crdt")
tab3m1q3 <- DifnDif3bis (a="welfare", b="consum.crdt")
tab3m1q4 <- DifnDif3bis (a="welfare", b="mrtg.crdt")
tab3m1q5 <- DifnDif3bis (a="unempl", b="consum.crdt")
tab3m1q6 <- DifnDif3bis (a="unempl", b="mrtg.crdt")

# Table 3, model 2 (high income)
tab3m2q1 <- DifnDif3bis_sub (a="tax", b="consum.crdt", subsetVar="income.cat"
                 , subsetValue="High Income")
tab3m2q2 <- DifnDif3bis_sub (a="tax", b="mrtg.crdt", subsetVar="income.cat"
                 , subsetValue="High Income")
tab3m2q3 <- DifnDif3bis_sub (a="welfare", b="consum.crdt", subsetVar="income.cat"
                 , subsetValue="High Income")
tab3m2q4 <- DifnDif3bis_sub (a="welfare", b="mrtg.crdt", subsetVar="income.cat"
                 , subsetValue="High Income")
tab3m2q5 <- DifnDif3bis_sub (a="unempl", b="consum.crdt", subsetVar="income.cat"
                 , subsetValue="High Income")
tab3m2q6 <- DifnDif3bis_sub (a="unempl", b="mrtg.crdt", subsetVar="income.cat"
                 , subsetValue="High Income")

# Table 3, model 3 (low income)
tab3m3q1 <- DifnDif3bis_sub (a="tax", b="consum.crdt", subsetVar="income.cat"
                 , subsetValue="Low Income")
tab3m3q2 <- DifnDif3bis_sub (a="tax", b="mrtg.crdt", subsetVar="income.cat"
                 , subsetValue="Low Income")
tab3m3q3 <- DifnDif3bis_sub (a="welfare", b="consum.crdt", subsetVar="income.cat"
                 , subsetValue="Low Income")
tab3m3q4 <- DifnDif3bis_sub (a="welfare", b="mrtg.crdt", subsetVar="income.cat"
                 , subsetValue="Low Income")
tab3m3q5 <- DifnDif3bis_sub (a="unempl", b="consum.crdt", subsetVar="income.cat"
                 , subsetValue="Low Income")
tab3m3q6 <- DifnDif3bis_sub (a="unempl", b="mrtg.crdt", subsetVar="income.cat"
                 , subsetValue="Low Income")

# Table 3, model 4 (low risk)
tab3m4q1 <- DifnDif3bis_sub (a="tax", b="consum.crdt", subsetVar="highUnemp"
                 , subsetValue="Low")
tab3m4q2 <- DifnDif3bis_sub (a="tax", b="mrtg.crdt", subsetVar="highUnemp"
                 , subsetValue="Low")
tab3m4q3 <- DifnDif3bis_sub (a="welfare", b="consum.crdt", subsetVar="highUnemp"
                 , subsetValue="Low")
tab3m4q4 <- DifnDif3bis_sub (a="welfare", b="mrtg.crdt", subsetVar="highUnemp"
                 , subsetValue="Low")
tab3m4q5 <- DifnDif3bis_sub (a="unempl", b="consum.crdt", subsetVar="highUnemp"
                 , subsetValue="Low")
tab3m4q6 <- DifnDif3bis_sub (a="unempl", b="mrtg.crdt", subsetVar="highUnemp"
                 , subsetValue="Low")

# Table 3, model 5 (high risk)
tab3m5q1 <- DifnDif3bis_sub (a="tax", b="consum.crdt", subsetVar="highUnemp"
                 , subsetValue="High")
tab3m5q2 <- DifnDif3bis_sub (a="tax", b="mrtg.crdt", subsetVar="highUnemp"
                 , subsetValue="High")
tab3m5q3 <- DifnDif3bis_sub (a="welfare", b="consum.crdt", subsetVar="highUnemp"
                 , subsetValue="High")
tab3m5q4 <- DifnDif3bis_sub (a="welfare", b="mrtg.crdt", subsetVar="highUnemp"
                 , subsetValue="High")
tab3m5q5 <- DifnDif3bis_sub (a="unempl", b="consum.crdt", subsetVar="highUnemp"
                 , subsetValue="High")
tab3m5q6 <- DifnDif3bis_sub (a="unempl", b="mrtg.crdt", subsetVar="highUnemp"
                 , subsetValue="High")

# Table 3
tab3mod1.cf <-  c(tab3m1q1$Coefficient, 
                  tab3m1q2$Coefficient, 
                  tab3m1q3$Coefficient, 
                  tab3m1q4$Coefficient, 
                  tab3m1q5$Coefficient, 
                  tab3m1q6$Coefficient)
tab3mod1.SE <-  c(tab3m1q1$SE, 
                  tab3m1q2$SE, 
                  tab3m1q3$SE, 
                  tab3m1q4$SE, 
                  tab3m1q5$SE, 
                  tab3m1q6$SE)
tab3mod1.tstat <-  c(tab3m1q1$t_stat, 
                  tab3m1q2$t_stat, 
                  tab3m1q3$t_stat, 
                  tab3m1q4$t_stat, 
                  tab3m1q5$t_stat, 
                  tab3m1q6$t_stat)
tab3mod1.st <- ifelse (tab3mod1.tstat > 2.576, "***", 
                       ifelse (tab3mod1.tstat <= 2.576 & tab3mod1.tstat > 1.96, "**",
                               ifelse (tab3mod1.tstat <= 1.96 & tab3mod1.tstat > 1.645, "*","")))

tab3mod2.cf <-  c(tab3m2q1$Coefficient, 
                  tab3m2q2$Coefficient, 
                  tab3m2q3$Coefficient, 
                  tab3m2q4$Coefficient, 
                  tab3m2q5$Coefficient, 
                  tab3m2q6$Coefficient)
tab3mod2.SE <-  c(tab3m2q1$SE, 
                  tab3m2q2$SE, 
                  tab3m2q3$SE, 
                  tab3m2q4$SE, 
                  tab3m2q5$SE, 
                  tab3m2q6$SE)
tab3mod2.tstat <-  c(tab3m2q1$t_stat, 
                     tab3m2q2$t_stat, 
                     tab3m2q3$t_stat, 
                     tab3m2q4$t_stat, 
                     tab3m2q5$t_stat, 
                     tab3m2q6$t_stat)
tab3mod2.st <- ifelse (tab3mod2.tstat > 2.576, "***", 
                       ifelse (tab3mod2.tstat <= 2.576 & tab3mod2.tstat > 1.96, "**",
                               ifelse (tab3mod2.tstat <= 1.96 & tab3mod2.tstat > 1.645, "*","")))


tab3mod3.cf <-  c(tab3m3q1$Coefficient, 
                  tab3m3q2$Coefficient, 
                  tab3m3q3$Coefficient, 
                  tab3m3q4$Coefficient, 
                  tab3m3q5$Coefficient, 
                  tab3m3q6$Coefficient)
tab3mod3.SE <-  c(tab3m3q1$SE, 
                  tab3m3q2$SE, 
                  tab3m3q3$SE, 
                  tab3m3q4$SE, 
                  tab3m3q5$SE, 
                  tab3m3q6$SE)
tab3mod3.tstat <-  c(tab3m3q1$t_stat, 
                     tab3m3q2$t_stat, 
                     tab3m3q3$t_stat, 
                     tab3m3q4$t_stat, 
                     tab3m3q5$t_stat, 
                     tab3m3q6$t_stat)
tab3mod3.st <- ifelse (tab3mod3.tstat > 2.576, "***", 
                       ifelse (tab3mod3.tstat <= 2.576 & tab3mod3.tstat > 1.96, "**",
                               ifelse (tab3mod3.tstat <= 1.96 & tab3mod3.tstat > 1.645, "*","")))


tab3mod4.cf <-  c(tab3m4q1$Coefficient, 
                  tab3m4q2$Coefficient, 
                  tab3m4q3$Coefficient, 
                  tab3m4q4$Coefficient, 
                  tab3m4q5$Coefficient, 
                  tab3m4q6$Coefficient)
tab3mod4.SE <-  c(tab3m4q1$SE, 
                  tab3m4q2$SE, 
                  tab3m4q3$SE, 
                  tab3m4q4$SE, 
                  tab3m4q5$SE, 
                  tab3m4q6$SE)
tab3mod4.tstat <-  c(tab3m4q1$t_stat, 
                     tab3m4q2$t_stat, 
                     tab3m4q3$t_stat, 
                     tab3m4q4$t_stat, 
                     tab3m4q5$t_stat, 
                     tab3m4q6$t_stat)
tab3mod4.st <- ifelse (tab3mod4.tstat > 2.576, "***", 
                       ifelse (tab3mod4.tstat <= 2.576 & tab3mod4.tstat > 1.96, "**",
                               ifelse (tab3mod4.tstat <= 1.96 & tab3mod4.tstat > 1.645, "*","")))

tab3mod5.cf <-  c(tab3m5q1$Coefficient, 
                  tab3m5q2$Coefficient, 
                  tab3m5q3$Coefficient, 
                  tab3m5q4$Coefficient, 
                  tab3m5q5$Coefficient, 
                  tab3m5q6$Coefficient)
tab3mod5.SE <-  c(tab3m5q1$SE, 
                  tab3m5q2$SE, 
                  tab3m5q3$SE, 
                  tab3m5q4$SE, 
                  tab3m5q5$SE, 
                  tab3m5q6$SE)
tab3mod5.tstat <-  c(tab3m5q1$t_stat, 
                     tab3m5q2$t_stat, 
                     tab3m5q3$t_stat, 
                     tab3m5q4$t_stat, 
                     tab3m5q5$t_stat, 
                     tab3m5q6$t_stat)
tab3mod5.st <- ifelse (tab3mod5.tstat > 2.576, "***", 
                       ifelse (tab3mod5.tstat <= 2.576 & tab3mod5.tstat > 1.96, "**",
                               ifelse (tab3mod5.tstat <= 1.96 & tab3mod5.tstat > 1.645, "*","")))


Table3 <- cbind (tab3mod1.coef = round(tab3mod1.cf, 3), 
                 tab3mod1.SE = round(tab3mod1.SE, 3), 
                 tab3mod1.st,
                 tab3mod2.coef = round(tab3mod2.cf, 3), 
                 tab3mod2.SE = round(tab3mod2.SE, 3), 
                 tab3mod2.st,
                 tab3mod3.coef = round(tab3mod3.cf, 3), 
                 tab3mod3.SE = round(tab3mod3.SE, 3), 
                 tab3mod3.st,
                 tab3mod4.coef = round(tab3mod4.cf, 3), 
                 tab3mod4.SE = round(tab3mod4.SE, 3), 
                 tab3mod4.st,
                 tab3mod5.coef = round(tab3mod5.cf, 3), 
                 tab3mod5.SE = round(tab3mod5.SE, 3), 
                 tab3mod5.st)

# N of respondents
tab3m1n <- length(unique(surveyData$ResponseId))
tab3m2n <- length(unique(surveyData[surveyData$income.cat == "High Income",]$ResponseId))
tab3m3n <- length(unique(surveyData[surveyData$income.cat == "Low Income",]$ResponseId))
tab3m4n <- length(unique(surveyData[surveyData$highUnemp == "Low",]$ResponseId))
tab3m5n <- length(unique(surveyData[surveyData$highUnemp == "High",]$ResponseId))

# add to table
Table3 <- rbind(Table3, 
                c(rep(tab3m1n,3),
                  rep(tab3m2n,3),
                  rep(tab3m3n,3),
                  rep(tab3m4n,3),
                  rep(tab3m5n,3)))

# define row names
var_names <- c("Income Tax x Credit Card Loans",
               "Income Tax x Home Loans",
               "Social Assistance x Credit Card Loans",
               "Social Assistance x Home Loans",
               "Unemployment Insurance x Credit Card Loans",
               "Unemployment Insurance x Home Loans",
               "Number of respondents")

rownames (Table3) <- var_names
print (Table3)




# Figure F.1: Marginal mean support by income groups ========================
surveyData$income.cat <- as.factor(surveyData$income.cat)
mm_byIncome <- cj(surveyData, fig2, id = ~ResponseId, estimate = "mm"
                  , by = ~income.cat
                  , feature_labels = list(mrtg.crdt = "Housing Loan"
                                          , consum.crdt = "Credit Card Loan"
                                          , welfare = "Social Assistance"
                                          , unempl = "Unemployment Support"
                                          , tax = "Income Tax"))

# pdf(file = "graphPlots/baseline_incGroups.pdf",
#     width = 7,
#     height = 5)
plot(mm_byIncome, group = "BY", vline = .5)
# dev.off()

# Figure F.2: Marginal mean support by risk groups ==========================
surveyData$unempRisk.cat <- as.factor(surveyData$unempRisk.cat)
mm_byEmplRisk <- cj(surveyData, fig2, id = ~ResponseId, estimate = "mm"
                    , by = ~unempRisk.cat
                    , feature_labels = list(mrtg.crdt = "Housing Loan"
                                            , consum.crdt = "Credit Card Loan"
                                            , welfare = "Social Assistance"
                                            , unempl = "Unemployment Support"
                                            , tax = "Income Tax"))

# pdf(file = "graphPlots/baseline_riskGroups.pdf",
#     width = 7,
#     height = 5)
plot(mm_byEmplRisk, group = "BY", vline = .5)
# dev.off()

# Table F.1: ================================================================
## subgroups: experience with loans
### Model 1: Experience with credit
tabf1m1q1 <- DifnDif3bis_sub (a="tax", b="consum.crdt", subsetVar="expLoan.cat"
                 , subsetValue="Experience with Credit")
tabf1m1q2 <- DifnDif3bis_sub (a="tax", b="mrtg.crdt", subsetVar="expLoan.cat"
                 , subsetValue="Experience with Credit")
tabf1m1q3 <- DifnDif3bis_sub (a="welfare", b="consum.crdt", subsetVar="expLoan.cat"
                 , subsetValue="Experience with Credit")
tabf1m1q4 <- DifnDif3bis_sub (a="welfare", b="mrtg.crdt", subsetVar="expLoan.cat"
                 , subsetValue="Experience with Credit")
tabf1m1q5 <- DifnDif3bis_sub (a="unempl", b="consum.crdt", subsetVar="expLoan.cat"
                 , subsetValue="Experience with Credit")
tabf1m1q6 <- DifnDif3bis_sub (a="unempl", b="mrtg.crdt", subsetVar="expLoan.cat"
                 , subsetValue="Experience with Credit")

### Model 2: No experience with credit
tabf1m2q1 <- DifnDif3bis_sub (a="tax", b="consum.crdt", subsetVar="expLoan.cat"
                 , subsetValue="No experience with Credit") 
tabf1m2q2 <- DifnDif3bis_sub (a="tax", b="mrtg.crdt", subsetVar="expLoan.cat"
                 , subsetValue="No experience with Credit") 
tabf1m2q3 <- DifnDif3bis_sub (a="welfare", b="consum.crdt", subsetVar="expLoan.cat"
                 , subsetValue="No experience with Credit") 
tabf1m2q4 <- DifnDif3bis_sub (a="welfare", b="mrtg.crdt", subsetVar="expLoan.cat"
                 , subsetValue="No experience with Credit") 
tabf1m2q5 <- DifnDif3bis_sub (a="unempl", b="consum.crdt", subsetVar="expLoan.cat"
                 , subsetValue="No experience with Credit") 
tabf1m2q6 <- DifnDif3bis_sub (a="unempl", b="mrtg.crdt", subsetVar="expLoan.cat"
                 , subsetValue="No experience with Credit") 

## subgroups: losses in financial crisis
### Model 3: Losses
tabf1m3q1 <- DifnDif3bis_sub (a="tax", b="consum.crdt", subsetVar="expFinCrisis.cat"
                 , subsetValue="Losses in financial crisis")
tabf1m3q2 <- DifnDif3bis_sub (a="tax", b="mrtg.crdt", subsetVar="expFinCrisis.cat"
                 , subsetValue="Losses in financial crisis")
tabf1m3q3 <- DifnDif3bis_sub (a="welfare", b="consum.crdt", subsetVar="expFinCrisis.cat"
                 , subsetValue="Losses in financial crisis")
tabf1m3q4 <- DifnDif3bis_sub (a="welfare", b="mrtg.crdt", subsetVar="expFinCrisis.cat"
                 , subsetValue="Losses in financial crisis")
tabf1m3q5 <- DifnDif3bis_sub (a="unempl", b="consum.crdt", subsetVar="expFinCrisis.cat"
                 , subsetValue="Losses in financial crisis")
tabf1m3q6 <- DifnDif3bis_sub (a="unempl", b="mrtg.crdt", subsetVar="expFinCrisis.cat"
                 , subsetValue="Losses in financial crisis")

### Model 4: No losses
tabf1m4q1 <- DifnDif3bis_sub (a="tax", b="consum.crdt", subsetVar="expFinCrisis.cat"
                 , subsetValue="No losses in financial crisis")
tabf1m4q2 <- DifnDif3bis_sub (a="tax", b="mrtg.crdt", subsetVar="expFinCrisis.cat"
                 , subsetValue="No losses in financial crisis")
tabf1m4q3 <- DifnDif3bis_sub (a="welfare", b="consum.crdt", subsetVar="expFinCrisis.cat"
                 , subsetValue="No losses in financial crisis")
tabf1m4q4 <- DifnDif3bis_sub (a="welfare", b="mrtg.crdt", subsetVar="expFinCrisis.cat"
                 , subsetValue="No losses in financial crisis")
tabf1m4q5 <- DifnDif3bis_sub (a="unempl", b="consum.crdt", subsetVar="expFinCrisis.cat"
                 , subsetValue="No losses in financial crisis")
tabf1m4q6 <- DifnDif3bis_sub (a="unempl", b="mrtg.crdt", subsetVar="expFinCrisis.cat"
                 , subsetValue="No losses in financial crisis")

tabf1mod1.cf <-  c(tabf1m1q1$Coefficient, 
                   tabf1m1q2$Coefficient, 
                   tabf1m1q3$Coefficient, 
                   tabf1m1q4$Coefficient, 
                   tabf1m1q5$Coefficient, 
                   tabf1m1q6$Coefficient)
tabf1mod1.SE <-  c(tabf1m1q1$SE, 
                  tabf1m1q2$SE, 
                  tabf1m1q3$SE, 
                  tabf1m1q4$SE, 
                  tabf1m1q5$SE, 
                  tabf1m1q6$SE)
tabf1mod1.tstat <-  c(tabf1m1q1$t_stat, 
                      tabf1m1q2$t_stat, 
                      tabf1m1q3$t_stat, 
                      tabf1m1q4$t_stat, 
                      tabf1m1q5$t_stat, 
                      tabf1m1q6$t_stat)
tabf1mod1.st <- ifelse (tabf1mod1.tstat > 2.576, "***", 
                       ifelse (tabf1mod1.tstat <= 2.576 & tabf1mod1.tstat > 1.96, "**",
                               ifelse (tabf1mod1.tstat <= 1.96 & tabf1mod1.tstat > 1.645, "*","")))

tabf1mod2.cf <-  c(tabf1m2q1$Coefficient, 
                   tabf1m2q2$Coefficient, 
                   tabf1m2q3$Coefficient, 
                   tabf1m2q4$Coefficient, 
                   tabf1m2q5$Coefficient, 
                   tabf1m2q6$Coefficient)
tabf1mod2.SE <-  c(tabf1m2q1$SE, 
                   tabf1m2q2$SE, 
                   tabf1m2q3$SE, 
                   tabf1m2q4$SE, 
                   tabf1m2q5$SE, 
                   tabf1m2q6$SE)
tabf1mod2.tstat <-  c(tabf1m2q1$t_stat, 
                      tabf1m2q2$t_stat, 
                      tabf1m2q3$t_stat, 
                      tabf1m2q4$t_stat, 
                      tabf1m2q5$t_stat, 
                      tabf1m2q6$t_stat)
tabf1mod2.st <- ifelse (tabf1mod2.tstat > 2.576, "***", 
                        ifelse (tabf1mod2.tstat <= 2.576 & tabf1mod2.tstat > 1.96, "**",
                                ifelse (tabf1mod2.tstat <= 1.96 & tabf1mod2.tstat > 1.645, "*","")))

tabf1mod3.cf <-  c(tabf1m3q1$Coefficient, 
                   tabf1m3q2$Coefficient, 
                   tabf1m3q3$Coefficient, 
                   tabf1m3q4$Coefficient, 
                   tabf1m3q5$Coefficient, 
                   tabf1m3q6$Coefficient)
tabf1mod3.SE <-  c(tabf1m3q1$SE, 
                   tabf1m3q2$SE, 
                   tabf1m3q3$SE, 
                   tabf1m3q4$SE, 
                   tabf1m3q5$SE, 
                   tabf1m3q6$SE)
tabf1mod3.tstat <-  c(tabf1m3q1$t_stat, 
                      tabf1m3q2$t_stat, 
                      tabf1m3q3$t_stat, 
                      tabf1m3q4$t_stat, 
                      tabf1m3q5$t_stat, 
                      tabf1m3q6$t_stat)
tabf1mod3.st <- ifelse (tabf1mod3.tstat > 2.576, "***", 
                        ifelse (tabf1mod3.tstat <= 2.576 & tabf1mod3.tstat > 1.96, "**",
                                ifelse (tabf1mod3.tstat <= 1.96 & tabf1mod3.tstat > 1.645, "*","")))


tabf1mod4.cf <-  c(tabf1m4q1$Coefficient, 
                   tabf1m4q2$Coefficient, 
                   tabf1m4q3$Coefficient, 
                   tabf1m4q4$Coefficient, 
                   tabf1m4q5$Coefficient, 
                   tabf1m4q6$Coefficient)
tabf1mod4.SE <-  c(tabf1m4q1$SE, 
                   tabf1m4q2$SE, 
                   tabf1m4q3$SE, 
                   tabf1m4q4$SE, 
                   tabf1m4q5$SE, 
                   tabf1m4q6$SE)
tabf1mod4.tstat <-  c(tabf1m4q1$t_stat, 
                      tabf1m4q2$t_stat, 
                      tabf1m4q3$t_stat, 
                      tabf1m4q4$t_stat, 
                      tabf1m4q5$t_stat, 
                      tabf1m4q6$t_stat)
tabf1mod4.st <- ifelse (tabf1mod4.tstat > 2.576, "***", 
                        ifelse (tabf1mod4.tstat <= 2.576 & tabf1mod4.tstat > 1.96, "**",
                                ifelse (tabf1mod4.tstat <= 1.96 & tabf1mod4.tstat > 1.645, "*","")))


TableF1 <- cbind (mod1.coef = round(tabf1mod1.cf, 3), 
                  mod1.SE = round(tabf1mod1.SE, 3), 
                 mod1.st = tabf1mod1.st,
                 mod2.coef = round(tabf1mod2.cf, 3), 
                 mod2.SE = round(tabf1mod2.SE, 3), 
                 mod2.st = tabf1mod2.st,
                 mod3.coef = round(tabf1mod3.cf, 3), 
                 mod3.SE = round(tabf1mod3.SE, 3), 
                 mod3.st = tabf1mod3.st,
                 mod4.coef = round(tabf1mod4.cf, 3), 
                 mod4.SE = round(tabf1mod4.SE, 3), 
                 mod4.st = tabf1mod4.st)

# N of respondents
tabf1m1n <- length(unique(surveyData[surveyData$expLoan.cat =="Experience with Credit",]$ResponseId))
tabf1m2n <- length(unique(surveyData[surveyData$expLoan.cat =="No experience with Credit",]$ResponseId))
tabf1m3n <- length(unique(surveyData[surveyData$expFinCrisis.cat == "Losses in financial crisis",]$ResponseId))
tabf1m4n <- length(unique(surveyData[surveyData$expFinCrisis.cat == "No losses in financial crisis",]$ResponseId))

# add to table
TableF1 <- rbind(TableF1, 
                c(rep(tabf1m1n,3),
                  rep(tabf1m2n,3),
                  rep(tabf1m3n,3),
                  rep(tabf1m4n,3)))

rownames (TableF1) <- var_names

print (TableF1)

# Table F.2: DID by Risk-Income Combination =================================

## subgroups: risk & income combined
### Model 1: high income-low risk
tabf2m1q1 <- DifnDif3bis_sub (a="tax", b="consum.crdt", subsetVar="incomeRisk.cat"
                 , subsetValue="High Income-Low Risk")
tabf2m1q2 <- DifnDif3bis_sub (a="tax", b="mrtg.crdt", subsetVar="incomeRisk.cat"
                , subsetValue="High Income-Low Risk") 
tabf2m1q3 <- DifnDif3bis_sub (a="welfare", b="consum.crdt", subsetVar="incomeRisk.cat"
                 , subsetValue="High Income-Low Risk")
tabf2m1q4 <- DifnDif3bis_sub (a="welfare", b="mrtg.crdt", subsetVar="incomeRisk.cat"
                 , subsetValue="High Income-Low Risk")
tabf2m1q5 <- DifnDif3bis_sub (a="unempl", b="consum.crdt", subsetVar="incomeRisk.cat"
                 , subsetValue="High Income-Low Risk")
tabf2m1q6 <- DifnDif3bis_sub (a="unempl", b="mrtg.crdt", subsetVar="incomeRisk.cat"
                 , subsetValue="High Income-Low Risk")

### Model 2: low income-high risk
tabf2m2q1 <- DifnDif3bis_sub (a="tax", b="consum.crdt", subsetVar="incomeRisk.cat"
                 , subsetValue="Low Income-High Risk")
tabf2m2q2 <- DifnDif3bis_sub (a="tax", b="mrtg.crdt", subsetVar="incomeRisk.cat"
                 , subsetValue="Low Income-High Risk")
tabf2m2q3 <- DifnDif3bis_sub (a="welfare", b="consum.crdt", subsetVar="incomeRisk.cat"
                 , subsetValue="Low Income-High Risk")
tabf2m2q4 <- DifnDif3bis_sub (a="welfare", b="mrtg.crdt", subsetVar="incomeRisk.cat"
                 , subsetValue="Low Income-High Risk")
tabf2m2q5 <- DifnDif3bis_sub (a="unempl", b="consum.crdt", subsetVar="incomeRisk.cat"
                 , subsetValue="Low Income-High Risk")
tabf2m2q6 <- DifnDif3bis_sub (a="unempl", b="mrtg.crdt", subsetVar="incomeRisk.cat"
                 , subsetValue="Low Income-High Risk")


tabf2mod1.cf <-  c(tabf2m1q1$Coefficient, 
                   tabf2m1q2$Coefficient, 
                   tabf2m1q3$Coefficient, 
                   tabf2m1q4$Coefficient, 
                   tabf2m1q5$Coefficient, 
                   tabf2m1q6$Coefficient)
tabf2mod1.SE <-  c(tabf2m1q1$SE, 
                   tabf2m1q2$SE, 
                   tabf2m1q3$SE, 
                   tabf2m1q4$SE, 
                   tabf2m1q5$SE, 
                   tabf2m1q6$SE)
tabf2mod1.tstat <-  c(tabf2m1q1$t_stat, 
                      tabf2m1q2$t_stat, 
                      tabf2m1q3$t_stat, 
                      tabf2m1q4$t_stat, 
                      tabf2m1q5$t_stat, 
                      tabf2m1q6$t_stat)
tabf2mod1.st <- ifelse (abs(tabf2mod1.tstat) > 2.576, "***", 
                        ifelse (abs(tabf2mod1.tstat) <= 2.576 & abs(tabf2mod1.tstat) > 1.96, "**",
                                ifelse (abs(tabf2mod1.tstat) <= 1.96 & abs(tabf2mod1.tstat) > 1.645, "*","")))

tabf2mod2.cf <-  c(tabf2m2q1$Coefficient, 
                   tabf2m2q2$Coefficient, 
                   tabf2m2q3$Coefficient, 
                   tabf2m2q4$Coefficient, 
                   tabf2m2q5$Coefficient, 
                   tabf2m2q6$Coefficient)
tabf2mod2.SE <-  c(tabf2m2q1$SE, 
                   tabf2m2q2$SE, 
                   tabf2m2q3$SE, 
                   tabf2m2q4$SE, 
                   tabf2m2q5$SE, 
                   tabf2m2q6$SE)
tabf2mod2.tstat <-  c(tabf2m2q1$t_stat, 
                      tabf2m2q2$t_stat, 
                      tabf2m2q3$t_stat, 
                      tabf2m2q4$t_stat, 
                      tabf2m2q5$t_stat, 
                      tabf2m2q6$t_stat)
tabf2mod2.st <- ifelse (abs(tabf2mod2.tstat) > 2.576, "***", 
                        ifelse (abs(tabf2mod2.tstat) <= 2.576 & abs(tabf2mod2.tstat) > 1.96, "**",
                                ifelse (abs(tabf2mod2.tstat) <= 1.96 & abs(tabf2mod2.tstat) > 1.645, "*","")))

TableF2 <- cbind (mod1.coef = round(tabf2mod1.cf, 3), 
                  mod1.SE = round(tabf2mod1.SE, 3), 
                  mod1.st = tabf2mod1.st,
                  mod2.coef = round(tabf2mod2.cf, 3), 
                  mod2.SE = round(tabf2mod2.SE, 3), 
                  mod2.st = tabf2mod2.st)

# N of respondents
tabf2m1n <- length(unique(surveyData[surveyData$incomeRisk.cat =="High Income-Low Risk",]$ResponseId))
tabf2m2n <- length(unique(surveyData[surveyData$incomeRisk.cat =="Low Income-High Risk",]$ResponseId))

# add to table
TableF2 <- rbind(TableF2, 
                 c(rep(tabf2m1n,3),
                   rep(tabf2m2n,3)))

rownames (TableF2) <- var_names

print (TableF2)

print (TableF2)


# Figure G.1: AMCE by typical/atyptical profile =============================

surveyData$atypicalProfile <- as.factor(surveyData$atypicalProfile)
amce_atypical <- cj(surveyData, fig2, id = ~ResponseId, estimate = "amce"
                    , by = ~atypicalProfile
                  , feature_labels = list(mrtg.crdt = "Housing Loan"
                                          , consum.crdt = "Credit Card Loan"
                                          , welfare = "Social Assistance"
                                          , unempl = "Unemployment Support"
                                          , tax = "Income Tax"))

# pdf(file = "graphPlots/AMCEtypical.pdf",
#     width = 7,
#     height = 5)
plot(amce_atypical, group = "BY", vline = 0)
# dev.off()

# Table G.1: DID excl atyptical profiles ====================================

# create new dataset for analysis
backupDta <- surveyData
surveyData <- surveyData[surveyData$atypicalProfile == F,]

# model 1 (full sample)
tabg1m1q1 <- DifnDif3bis (a="tax", b="consum.crdt")
tabg1m1q2 <- DifnDif3bis (a="tax", b="mrtg.crdt")
tabg1m1q3 <- DifnDif3bis (a="welfare", b="consum.crdt")
tabg1m1q4 <- DifnDif3bis (a="welfare", b="mrtg.crdt")
tabg1m1q5 <- DifnDif3bis (a="unempl", b="consum.crdt")
tabg1m1q6 <- DifnDif3bis (a="unempl", b="mrtg.crdt")

# model 2 (high income)
tabg1m2q1 <- DifnDif3bis_sub (a="tax", b="consum.crdt", subsetVar="income.cat"
                             , subsetValue="High Income")
tabg1m2q2 <- DifnDif3bis_sub (a="tax", b="mrtg.crdt", subsetVar="income.cat"
                             , subsetValue="High Income")
tabg1m2q3 <- DifnDif3bis_sub (a="welfare", b="consum.crdt", subsetVar="income.cat"
                             , subsetValue="High Income")
tabg1m2q4 <- DifnDif3bis_sub (a="welfare", b="mrtg.crdt", subsetVar="income.cat"
                             , subsetValue="High Income")
tabg1m2q5 <- DifnDif3bis_sub (a="unempl", b="consum.crdt", subsetVar="income.cat"
                             , subsetValue="High Income")
tabg1m2q6 <- DifnDif3bis_sub (a="unempl", b="mrtg.crdt", subsetVar="income.cat"
                             , subsetValue="High Income")

# model 3 (low income)
tabg1m3q1 <- DifnDif3bis_sub (a="tax", b="consum.crdt", subsetVar="income.cat"
                             , subsetValue="Low Income")
tabg1m3q2 <- DifnDif3bis_sub (a="tax", b="mrtg.crdt", subsetVar="income.cat"
                             , subsetValue="Low Income")
tabg1m3q3 <- DifnDif3bis_sub (a="welfare", b="consum.crdt", subsetVar="income.cat"
                             , subsetValue="Low Income")
tabg1m3q4 <- DifnDif3bis_sub (a="welfare", b="mrtg.crdt", subsetVar="income.cat"
                             , subsetValue="Low Income")
tabg1m3q5 <- DifnDif3bis_sub (a="unempl", b="consum.crdt", subsetVar="income.cat"
                             , subsetValue="Low Income")
tabg1m3q6 <- DifnDif3bis_sub (a="unempl", b="mrtg.crdt", subsetVar="income.cat"
                             , subsetValue="Low Income")

# model 4 (low risk)
tabg1m4q1 <- DifnDif3bis_sub (a="tax", b="consum.crdt", subsetVar="highUnemp"
                             , subsetValue="Low")
tabg1m4q2 <- DifnDif3bis_sub (a="tax", b="mrtg.crdt", subsetVar="highUnemp"
                             , subsetValue="Low")
tabg1m4q3 <- DifnDif3bis_sub (a="welfare", b="consum.crdt", subsetVar="highUnemp"
                             , subsetValue="Low")
tabg1m4q4 <- DifnDif3bis_sub (a="welfare", b="mrtg.crdt", subsetVar="highUnemp"
                             , subsetValue="Low")
tabg1m4q5 <- DifnDif3bis_sub (a="unempl", b="consum.crdt", subsetVar="highUnemp"
                             , subsetValue="Low")
tabg1m4q6 <- DifnDif3bis_sub (a="unempl", b="mrtg.crdt", subsetVar="highUnemp"
                             , subsetValue="Low")

# model 5 (high risk)
tabg1m5q1 <- DifnDif3bis_sub (a="tax", b="consum.crdt", subsetVar="highUnemp"
                             , subsetValue="High")
tabg1m5q2 <- DifnDif3bis_sub (a="tax", b="mrtg.crdt", subsetVar="highUnemp"
                             , subsetValue="High")
tabg1m5q3 <- DifnDif3bis_sub (a="welfare", b="consum.crdt", subsetVar="highUnemp"
                             , subsetValue="High")
tabg1m5q4 <- DifnDif3bis_sub (a="welfare", b="mrtg.crdt", subsetVar="highUnemp"
                             , subsetValue="High")
tabg1m5q5 <- DifnDif3bis_sub (a="unempl", b="consum.crdt", subsetVar="highUnemp"
                             , subsetValue="High")
tabg1m5q6 <- DifnDif3bis_sub (a="unempl", b="mrtg.crdt", subsetVar="highUnemp"
                             , subsetValue="High")


tabg1mod1.cf <-  c(tabg1m1q1$Coefficient, 
                  tabg1m1q2$Coefficient, 
                  tabg1m1q3$Coefficient, 
                  tabg1m1q4$Coefficient, 
                  tabg1m1q5$Coefficient, 
                  tabg1m1q6$Coefficient)
tabg1mod1.SE <-  c(tabg1m1q1$SE, 
                  tabg1m1q2$SE, 
                  tabg1m1q3$SE, 
                  tabg1m1q4$SE, 
                  tabg1m1q5$SE, 
                  tabg1m1q6$SE)
tabg1mod1.tstat <-  c(tabg1m1q1$t_stat, 
                     tabg1m1q2$t_stat, 
                     tabg1m1q3$t_stat, 
                     tabg1m1q4$t_stat, 
                     tabg1m1q5$t_stat, 
                     tabg1m1q6$t_stat)
tabg1mod1.st <- ifelse (tabg1mod1.tstat > 2.576, "***", 
                       ifelse (tabg1mod1.tstat <= 2.576 & tabg1mod1.tstat > 1.96, "**",
                               ifelse (tabg1mod1.tstat <= 1.96 & tabg1mod1.tstat > 1.645, "*","")))

tabg1mod2.cf <-  c(tabg1m2q1$Coefficient, 
                  tabg1m2q2$Coefficient, 
                  tabg1m2q3$Coefficient, 
                  tabg1m2q4$Coefficient, 
                  tabg1m2q5$Coefficient, 
                  tabg1m2q6$Coefficient)
tabg1mod2.SE <-  c(tabg1m2q1$SE, 
                  tabg1m2q2$SE, 
                  tabg1m2q3$SE, 
                  tabg1m2q4$SE, 
                  tabg1m2q5$SE, 
                  tabg1m2q6$SE)
tabg1mod2.tstat <-  c(tabg1m2q1$t_stat, 
                     tabg1m2q2$t_stat, 
                     tabg1m2q3$t_stat, 
                     tabg1m2q4$t_stat, 
                     tabg1m2q5$t_stat, 
                     tabg1m2q6$t_stat)
tabg1mod2.st <- ifelse (tabg1mod2.tstat > 2.576, "***", 
                       ifelse (tabg1mod2.tstat <= 2.576 & tabg1mod2.tstat > 1.96, "**",
                               ifelse (tabg1mod2.tstat <= 1.96 & tabg1mod2.tstat > 1.645, "*","")))


tabg1mod3.cf <-  c(tabg1m3q1$Coefficient, 
                  tabg1m3q2$Coefficient, 
                  tabg1m3q3$Coefficient, 
                  tabg1m3q4$Coefficient, 
                  tabg1m3q5$Coefficient, 
                  tabg1m3q6$Coefficient)
tabg1mod3.SE <-  c(tabg1m3q1$SE, 
                  tabg1m3q2$SE, 
                  tabg1m3q3$SE, 
                  tabg1m3q4$SE, 
                  tabg1m3q5$SE, 
                  tabg1m3q6$SE)
tabg1mod3.tstat <-  c(tabg1m3q1$t_stat, 
                     tabg1m3q2$t_stat, 
                     tabg1m3q3$t_stat, 
                     tabg1m3q4$t_stat, 
                     tabg1m3q5$t_stat, 
                     tabg1m3q6$t_stat)
tabg1mod3.st <- ifelse (tabg1mod3.tstat > 2.576, "***", 
                       ifelse (tabg1mod3.tstat <= 2.576 & tabg1mod3.tstat > 1.96, "**",
                               ifelse (tabg1mod3.tstat <= 1.96 & tabg1mod3.tstat > 1.645, "*","")))


tabg1mod4.cf <-  c(tabg1m4q1$Coefficient, 
                  tabg1m4q2$Coefficient, 
                  tabg1m4q3$Coefficient, 
                  tabg1m4q4$Coefficient, 
                  tabg1m4q5$Coefficient, 
                  tabg1m4q6$Coefficient)
tabg1mod4.SE <-  c(tabg1m4q1$SE, 
                  tabg1m4q2$SE, 
                  tabg1m4q3$SE, 
                  tabg1m4q4$SE, 
                  tabg1m4q5$SE, 
                  tabg1m4q6$SE)
tabg1mod4.tstat <-  c(tabg1m4q1$t_stat, 
                     tabg1m4q2$t_stat, 
                     tabg1m4q3$t_stat, 
                     tabg1m4q4$t_stat, 
                     tabg1m4q5$t_stat, 
                     tabg1m4q6$t_stat)
tabg1mod4.st <- ifelse (tabg1mod4.tstat > 2.576, "***", 
                       ifelse (tabg1mod4.tstat <= 2.576 & tabg1mod4.tstat > 1.96, "**",
                               ifelse (tabg1mod4.tstat <= 1.96 & tabg1mod4.tstat > 1.645, "*","")))

tabg1mod5.cf <-  c(tabg1m5q1$Coefficient, 
                  tabg1m5q2$Coefficient, 
                  tabg1m5q3$Coefficient, 
                  tabg1m5q4$Coefficient, 
                  tabg1m5q5$Coefficient, 
                  tabg1m5q6$Coefficient)
tabg1mod5.SE <-  c(tabg1m5q1$SE, 
                  tabg1m5q2$SE, 
                  tabg1m5q3$SE, 
                  tabg1m5q4$SE, 
                  tabg1m5q5$SE, 
                  tabg1m5q6$SE)
tabg1mod5.tstat <-  c(tabg1m5q1$t_stat, 
                     tabg1m5q2$t_stat, 
                     tabg1m5q3$t_stat, 
                     tabg1m5q4$t_stat, 
                     tabg1m5q5$t_stat, 
                     tabg1m5q6$t_stat)
tabg1mod5.st <- ifelse (tabg1mod5.tstat > 2.576, "***", 
                       ifelse (tabg1mod5.tstat <= 2.576 & tabg1mod5.tstat > 1.96, "**",
                               ifelse (tabg1mod5.tstat <= 1.96 & tabg1mod5.tstat > 1.645, "*","")))


TableG1 <- cbind (mod1.coef = round(tabg1mod1.cf, 3), 
                 mod1.SE = round(tabg1mod1.SE, 3), 
                 mod5.st = tabg1mod1.st,
                 mod2.coef = round(tabg1mod2.cf, 3), 
                 mod2.SE = round(tabg1mod2.SE, 3), 
                 mod5.st = tabg1mod2.st,
                 mod3.coef = round(tabg1mod3.cf, 3), 
                 mod3.SE = round(tabg1mod3.SE, 3), 
                 mod5.st = tabg1mod3.st,
                 mod4.coef = round(tabg1mod4.cf, 3), 
                 mod4.SE = round(tabg1mod4.SE, 3), 
                 mod5.st = tabg1mod4.st,
                 mod5.coef = round(tabg1mod5.cf, 3), 
                 mod5.SE = round(tabg1mod5.SE, 3), 
                 mod5.st = tabg1mod5.st)

# N of respondents
tabg1m1n <- length(unique(surveyData$ResponseId))
tabg1m2n <- length(unique(surveyData[surveyData$income.cat == "High Income",]$ResponseId))
tabg1m3n <- length(unique(surveyData[surveyData$income.cat == "Low Income",]$ResponseId))
tabg1m4n <- length(unique(surveyData[surveyData$highUnemp == "Low",]$ResponseId))
tabg1m5n <- length(unique(surveyData[surveyData$highUnemp == "High",]$ResponseId))

# add to table
TableG1 <- rbind(TableG1, 
                c(rep(tabg1m1n,3),
                  rep(tabg1m2n,3),
                  rep(tabg1m3n,3),
                  rep(tabg1m4n,3),
                  rep(tabg1m5n,3)))

rownames(TableG1) <- var_names

print (TableG1)

# recreate old dataset
surveyData <- backupDta

# Table G.2: DID by ideology ================================================

# recreate old dataset
surveyData <- backupDta

### Model 1: left
tabg2m1q1 <- DifnDif3bis_sub (a="tax", b="consum.crdt", subsetVar="lr.cat"
                 , subsetValue="left")
tabg2m1q2 <- DifnDif3bis_sub (a="tax", b="mrtg.crdt", subsetVar="lr.cat"
                 , subsetValue="left")
tabg2m1q3 <- DifnDif3bis_sub (a="welfare", b="consum.crdt", subsetVar="lr.cat"
                 , subsetValue="left")
tabg2m1q4 <- DifnDif3bis_sub (a="welfare", b="mrtg.crdt", subsetVar="lr.cat"
                 , subsetValue="left")
tabg2m1q5 <- DifnDif3bis_sub (a="unempl", b="consum.crdt", subsetVar="lr.cat"
                 , subsetValue="left")
tabg2m1q6 <- DifnDif3bis_sub (a="unempl", b="mrtg.crdt", subsetVar="lr.cat"
                 , subsetValue="left")

### Model 2: right
tabg2m2q1 <- DifnDif3bis_sub (a="tax", b="consum.crdt", subsetVar="lr.cat"
                 , subsetValue="right")
tabg2m2q2 <- DifnDif3bis_sub (a="tax", b="mrtg.crdt", subsetVar="lr.cat"
                 , subsetValue="right")
tabg2m2q3 <- DifnDif3bis_sub (a="welfare", b="consum.crdt", subsetVar="lr.cat"
                 , subsetValue="right")
tabg2m2q4 <- DifnDif3bis_sub (a="welfare", b="mrtg.crdt", subsetVar="lr.cat"
                 , subsetValue="right")
tabg2m2q5 <- DifnDif3bis_sub (a="unempl", b="consum.crdt", subsetVar="lr.cat"
                 , subsetValue="right")
tabg2m2q6 <- DifnDif3bis_sub (a="unempl", b="mrtg.crdt", subsetVar="lr.cat"
                 , subsetValue="right")

tabg2mod1.cf <-  c(tabg2m1q1$Coefficient, 
                   tabg2m1q2$Coefficient, 
                   tabg2m1q3$Coefficient, 
                   tabg2m1q4$Coefficient, 
                   tabg2m1q5$Coefficient, 
                   tabg2m1q6$Coefficient)
tabg2mod1.SE <-  c(tabg2m1q1$SE, 
                   tabg2m1q2$SE, 
                   tabg2m1q3$SE, 
                   tabg2m1q4$SE, 
                   tabg2m1q5$SE, 
                   tabg2m1q6$SE)
tabg2mod1.tstat <-  c(tabg2m1q1$t_stat, 
                      tabg2m1q2$t_stat, 
                      tabg2m1q3$t_stat, 
                      tabg2m1q4$t_stat, 
                      tabg2m1q5$t_stat, 
                      tabg2m1q6$t_stat)
tabg2mod1.st <- ifelse (abs(tabg2mod1.tstat) > 2.576, "***", 
                        ifelse (abs(tabg2mod1.tstat) <= 2.576 & abs(tabg2mod1.tstat) > 1.96, "**",
                                ifelse (abs(tabg2mod1.tstat) <= 1.96 & abs(tabg2mod1.tstat) > 1.645, "*","")))

tabg2mod2.cf <-  c(tabg2m2q1$Coefficient, 
                   tabg2m2q2$Coefficient, 
                   tabg2m2q3$Coefficient, 
                   tabg2m2q4$Coefficient, 
                   tabg2m2q5$Coefficient, 
                   tabg2m2q6$Coefficient)
tabg2mod2.SE <-  c(tabg2m2q1$SE, 
                   tabg2m2q2$SE, 
                   tabg2m2q3$SE, 
                   tabg2m2q4$SE, 
                   tabg2m2q5$SE, 
                   tabg2m2q6$SE)
tabg2mod2.tstat <-  c(tabg2m2q1$t_stat, 
                      tabg2m2q2$t_stat, 
                      tabg2m2q3$t_stat, 
                      tabg2m2q4$t_stat, 
                      tabg2m2q5$t_stat, 
                      tabg2m2q6$t_stat)
tabg2mod2.st <- ifelse (abs(tabg2mod2.tstat) > 2.576, "***", 
                        ifelse (abs(tabg2mod2.tstat) <= 2.576 & abs(tabg2mod2.tstat) > 1.96, "**",
                                ifelse (abs(tabg2mod2.tstat) <= 1.96 & abs(tabg2mod2.tstat) > 1.645, "*","")))


TableG2 <- cbind (mod1.coef = round(tabg2mod1.cf, 3), 
                  mod1.SE = round(tabg2mod1.SE, 3), 
                  mod1.st = tabg2mod1.st,
                  mod2.coef = round(tabg2mod2.cf, 3), 
                  mod2.SE = round(tabg2mod2.SE, 3), 
                  mod2.st = tabg2mod2.st)

# N of respondents
tabg2m1n <- length(unique(surveyData[surveyData$lr.cat == "left",]$ResponseId))
tabg2m2n <- length(unique(surveyData[surveyData$lr.cat == "right",]$ResponseId))

# add to table
TableG2 <- rbind(TableG2, 
                 c(rep(tabg2m1n,3),
                   rep(tabg2m2n,3)))

rownames(TableG2) <- var_names

print (TableG2)

# Table G.3: DID by education ===============================================

### Model 1: higher education
tabg3m1q1 <- DifnDif3bis_sub (a="tax", b="consum.crdt", subsetVar="educ.cat"
                 , subsetValue="Higher Education")
tabg3m1q2 <- DifnDif3bis_sub (a="tax", b="mrtg.crdt", subsetVar="educ.cat"
                 , subsetValue="Higher Education")
tabg3m1q3 <- DifnDif3bis_sub (a="welfare", b="consum.crdt", subsetVar="educ.cat"
                 , subsetValue="Higher Education")
tabg3m1q4 <- DifnDif3bis_sub (a="welfare", b="mrtg.crdt", subsetVar="educ.cat"
                 , subsetValue="Higher Education")
tabg3m1q5 <- DifnDif3bis_sub (a="unempl", b="consum.crdt", subsetVar="educ.cat"
                 , subsetValue="Higher Education")
tabg3m1q6 <- DifnDif3bis_sub (a="unempl", b="mrtg.crdt", subsetVar="educ.cat"
                 , subsetValue="Higher Education")

### Model 2: no higher education
tabg3m2q1 <- DifnDif3bis_sub (a="tax", b="consum.crdt", subsetVar="educ.cat"
                 , subsetValue="No Higher Education") 
tabg3m2q2 <- DifnDif3bis_sub (a="tax", b="mrtg.crdt", subsetVar="educ.cat"
                , subsetValue="No Higher Education") 
tabg3m2q3 <- DifnDif3bis_sub (a="welfare", b="consum.crdt", subsetVar="educ.cat"
                 , subsetValue="No Higher Education")
tabg3m2q4 <- DifnDif3bis_sub (a="welfare", b="mrtg.crdt", subsetVar="educ.cat"
                 , subsetValue="No Higher Education")
tabg3m2q5 <- DifnDif3bis_sub (a="unempl", b="consum.crdt", subsetVar="educ.cat"
                 , subsetValue="No Higher Education")
tabg3m2q6 <- DifnDif3bis_sub (a="unempl", b="mrtg.crdt", subsetVar="educ.cat"
                 , subsetValue="No Higher Education")

tabg3mod1.cf <-  c(tabg3m1q1$Coefficient, 
                   tabg3m1q2$Coefficient, 
                   tabg3m1q3$Coefficient, 
                   tabg3m1q4$Coefficient, 
                   tabg3m1q5$Coefficient, 
                   tabg3m1q6$Coefficient)
tabg3mod1.SE <-  c(tabg3m1q1$SE, 
                   tabg3m1q2$SE, 
                   tabg3m1q3$SE, 
                   tabg3m1q4$SE, 
                   tabg3m1q5$SE, 
                   tabg3m1q6$SE)
tabg3mod1.tstat <-  c(tabg3m1q1$t_stat, 
                      tabg3m1q2$t_stat, 
                      tabg3m1q3$t_stat, 
                      tabg3m1q4$t_stat, 
                      tabg3m1q5$t_stat, 
                      tabg3m1q6$t_stat)
tabg3mod1.st <- ifelse (abs(tabg3mod1.tstat) > 2.576, "***", 
                        ifelse (abs(tabg3mod1.tstat) <= 2.576 & abs(tabg3mod1.tstat) > 1.96, "**",
                                ifelse (abs(tabg3mod1.tstat) <= 1.96 & abs(tabg3mod1.tstat) > 1.645, "*","")))

tabg3mod2.cf <-  c(tabg3m2q1$Coefficient, 
                   tabg3m2q2$Coefficient, 
                   tabg3m2q3$Coefficient, 
                   tabg3m2q4$Coefficient, 
                   tabg3m2q5$Coefficient, 
                   tabg3m2q6$Coefficient)
tabg3mod2.SE <-  c(tabg3m2q1$SE, 
                   tabg3m2q2$SE, 
                   tabg3m2q3$SE, 
                   tabg3m2q4$SE, 
                   tabg3m2q5$SE, 
                   tabg3m2q6$SE)
tabg3mod2.tstat <-  c(tabg3m2q1$t_stat, 
                      tabg3m2q2$t_stat, 
                      tabg3m2q3$t_stat, 
                      tabg3m2q4$t_stat, 
                      tabg3m2q5$t_stat, 
                      tabg3m2q6$t_stat)
tabg3mod2.st <- ifelse (abs(tabg3mod2.tstat) > 2.576, "***", 
                        ifelse (abs(tabg3mod2.tstat) <= 2.576 & abs(tabg3mod2.tstat) > 1.96, "**",
                                ifelse (abs(tabg3mod2.tstat) <= 1.96 & abs(tabg3mod2.tstat) > 1.645, "*","")))


TableG3 <- cbind (mod1.coef = round(tabg3mod1.cf, 3), 
                  mod1.SE = round(tabg3mod1.SE, 3), 
                  mod1.st = tabg3mod1.st,
                  mod2.coef = round(tabg3mod2.cf, 3), 
                  mod2.SE = round(tabg3mod2.SE, 3), 
                  mod2.st = tabg3mod2.st)

# N of respondents
tabg3m1n <- length(unique(surveyData[surveyData$educ.cat == "Higher Education",]$ResponseId))
tabg3m2n <- length(unique(surveyData[surveyData$educ.cat == "No Higher Education",]$ResponseId))

# add to table
TableG3 <- rbind(TableG3, 
                 c(rep(tabg3m1n,3),
                   rep(tabg3m2n,3)))

rownames(TableG3) <- var_names

print (TableG3)

# Figure G.2: AMCE by attention checks ======================================

surveyData$round.cat <- as.factor(surveyData$round.cat)
amce_byAttention <- cj(surveyData, fig2, id = ~ResponseId, estimate = "amce"
                       , by = ~round.cat
                       , feature_labels = list(mrtg.crdt = "Housing Loan"
                                               , consum.crdt = "Credit Card Loan"
                                               , welfare = "Social Assistance"
                                               , unempl = "Unemployment Support"
                                               , tax = "Income Tax"))
plot(amce_byAttention, group = "BY", vline = 0)

# Table G.4: DID by alternative risk measure ================================

### Model 1: low risk
tabg4m1q1 <- DifnDif3bis_sub (a="tax", b="consum.crdt", subsetVar="unempRisk.cat"
                 , subsetValue="Low Unemployment Risk") # line 1
tabg4m1q2 <- DifnDif3bis_sub (a="tax", b="mrtg.crdt", subsetVar="unempRisk.cat"
                 , subsetValue="Low Unemployment Risk") # line 2
tabg4m1q3 <- DifnDif3bis_sub (a="welfare", b="consum.crdt", subsetVar="unempRisk.cat"
                 , subsetValue="Low Unemployment Risk") # line 3
tabg4m1q4 <- DifnDif3bis_sub (a="welfare", b="mrtg.crdt", subsetVar="unempRisk.cat"
                 , subsetValue="Low Unemployment Risk") # line 4
tabg4m1q5 <- DifnDif3bis_sub (a="unempl", b="consum.crdt", subsetVar="unempRisk.cat"
                 , subsetValue="Low Unemployment Risk") # line 5
tabg4m1q6 <- DifnDif3bis_sub (a="unempl", b="mrtg.crdt", subsetVar="unempRisk.cat"
                 , subsetValue="Low Unemployment Risk") # line 6

### Model 2: high risk
tabg4m2q1 <- DifnDif3bis_sub (a="tax", b="consum.crdt", subsetVar="unempRisk.cat"
                 , subsetValue="High Unemployment Risk") # line 1
tabg4m2q2 <- DifnDif3bis_sub (a="tax", b="mrtg.crdt", subsetVar="unempRisk.cat"
                 , subsetValue="High Unemployment Risk") # line 2
tabg4m2q3 <- DifnDif3bis_sub (a="welfare", b="consum.crdt", subsetVar="unempRisk.cat"
                 , subsetValue="High Unemployment Risk") # line 3
tabg4m2q4 <- DifnDif3bis_sub (a="welfare", b="mrtg.crdt", subsetVar="unempRisk.cat"
                 , subsetValue="High Unemployment Risk") # line 4
tabg4m2q5 <- DifnDif3bis_sub (a="unempl", b="consum.crdt", subsetVar="unempRisk.cat"
                 , subsetValue="High Unemployment Risk") # line 5
tabg4m2q6 <- DifnDif3bis_sub (a="unempl", b="mrtg.crdt", subsetVar="unempRisk.cat"
                 , subsetValue="High Unemployment Risk") # line 6

tabg4mod1.cf <-  c(tabg4m1q1$Coefficient, 
                   tabg4m1q2$Coefficient, 
                   tabg4m1q3$Coefficient, 
                   tabg4m1q4$Coefficient, 
                   tabg4m1q5$Coefficient, 
                   tabg4m1q6$Coefficient)
tabg4mod1.SE <-  c(tabg4m1q1$SE, 
                   tabg4m1q2$SE, 
                   tabg4m1q3$SE, 
                   tabg4m1q4$SE, 
                   tabg4m1q5$SE, 
                   tabg4m1q6$SE)
tabg4mod1.tstat <-  c(tabg4m1q1$t_stat, 
                      tabg4m1q2$t_stat, 
                      tabg4m1q3$t_stat, 
                      tabg4m1q4$t_stat, 
                      tabg4m1q5$t_stat, 
                      tabg4m1q6$t_stat)
tabg4mod1.st <- ifelse (abs(tabg4mod1.tstat) > 2.576, "***", 
                        ifelse (abs(tabg4mod1.tstat) <= 2.576 & abs(tabg4mod1.tstat) > 1.96, "**",
                                ifelse (abs(tabg4mod1.tstat) <= 1.96 & abs(tabg4mod1.tstat) > 1.645, "*","")))

tabg4mod2.cf <-  c(tabg4m2q1$Coefficient, 
                   tabg4m2q2$Coefficient, 
                   tabg4m2q3$Coefficient, 
                   tabg4m2q4$Coefficient, 
                   tabg4m2q5$Coefficient, 
                   tabg4m2q6$Coefficient)
tabg4mod2.SE <-  c(tabg4m2q1$SE, 
                   tabg4m2q2$SE, 
                   tabg4m2q3$SE, 
                   tabg4m2q4$SE, 
                   tabg4m2q5$SE, 
                   tabg4m2q6$SE)
tabg4mod2.tstat <-  c(tabg4m2q1$t_stat, 
                      tabg4m2q2$t_stat, 
                      tabg4m2q3$t_stat, 
                      tabg4m2q4$t_stat, 
                      tabg4m2q5$t_stat, 
                      tabg4m2q6$t_stat)
tabg4mod2.st <- ifelse (abs(tabg4mod2.tstat) > 2.576, "***", 
                        ifelse (abs(tabg4mod2.tstat) <= 2.576 & abs(tabg4mod2.tstat) > 1.96, "**",
                                ifelse (abs(tabg4mod2.tstat) <= 1.96 & abs(tabg4mod2.tstat) > 1.645, "*","")))


TableG4 <- cbind (mod1.coef = round(tabg4mod1.cf, 3), 
                  mod1.SE = round(tabg4mod1.SE, 3), 
                  mod1.st = tabg4mod1.st,
                  mod2.coef = round(tabg4mod2.cf, 3), 
                  mod2.SE = round(tabg4mod2.SE, 3), 
                  mod2.st = tabg4mod2.st)

# N of respondents
tabg4m1n <- length(unique(surveyData[surveyData$unempRisk.cat == "Low Unemployment Risk",]$ResponseId))
tabg4m2n <- length(unique(surveyData[surveyData$unempRisk.cat == "High Unemployment Risk",]$ResponseId))

# add to table
TableG4 <- rbind(TableG4, 
                 c(rep(tabg4m1n,3),
                   rep(tabg4m2n,3)))

rownames(TableG4) <- var_names

print (TableG4)

